Semiparametric Paradigm Free Mapping: Automatic detection and characterization of fMRI BOLD responses and physiological fluctuations without prior information

Introduction: In recent work we showed that by means of sparse estimation techniques the spatial and temporal evolution of single-trial BOLD responses can be automatically detected without any prior knowledge of the stimulus timing and without thresholding: paradigm free mapping (PFM) [1]. However, fMRI time series also contain physiological and instrumental fluctuations which can hinder the detection of BOLD responses associated to neuronal activity [2]. Physiological fluctuations can be removed prior to PFM via high-pass filtering, or by RETROICOR, RVT or RVHRCOR [3,4,5], but these techniques must be employed in a preprocessing stage and require the additional recording of physiological respiratory and cardiac waveforms. Here, extending on our previous work, we present a novel technique which by decomposing the fMRI signal enables automatic detection of fMRI BOLD responses without prior stimulus information and automatic fitting of significant frequency fluctuations present in the signal, such as non-neuronal cardiac and respiratory fluctuations (semiparametric PFM, sPFM). This technique is based on a semiparametric linear representation of the fMRI signal [6,7] which is recursively fitted using a morphological component analysis algorithm [8]. The feasibility of this technique was evaluated in simulations and real fMRI data acquired at 7T, and its performance validated to RETROICOR. Semiparametric Parametric Free Mapping: Let us assume the fMRI signal at time point n can be written as yn = b(tn)+g(tn) +e(tn), where b(tn) represents the BOLD responses, g(tn) an arbitrary and unknown function and e(tn) a noise term. Assuming a linear model for the BOLD response, b(tn) can be described as the convolution of a stimulus-related signal s(t) with the BOLD response h(t). For N fMRI volumes, this component can be written in vectorial form as b=Hs, where H is a convolution matrix defined according to the BOLD response, and b and s are vectors of size N. Similarly, let us expand g as a linear combination of terms of the discrete cosine and sine (DCST) transforms, g = Φα, where Φ is a orthonomal matrix with the DCST terms and α is a vector with the amplitude of the components. Hence, the model signal can be written as y = Hs + Φα +e. Under this model, the vectors s and α can be estimated by minimizing a L1-regularization problem which minimizes the squared sum of the residuals and the L1-norms of both estimates. This optimization problem can be solved by morphological component analysis which alternates the estimation of one of the vectors (s and α) with the Basis Pursuit Denoising (BPDN) estimator whereas keeping constant the other estimate [8]. Methods: Simulations: fMRI time series were simulated by generating random stimulus waveforms (single-trial stimuli with random onset, ISI, and duration ranging from 2 to 6 s) which were convolved with an HRF with random parameters (mean=SPM Canonical HRF, std=0.2*mean). Subsequently, random Gaussian noise (CNR=2..4) and sinusoidal signals at 1Hz and 0.3 Hz (and harmonics) to simulate cardiac and respiratory fluctuations were added to the ideal waveforms. Different TR values were also simulated to study the effect of the Nyquist sampling criteria of physiological fluctuations. Real Data: Two fMRI datasets acquired on a 7T Philips scanner during a visual-motor paradigm were used (single-shot GE EPI, TR=2s, TE=30ms, 2mm isotropic resolution, 1 slice, 342 scans). Datasets were only corrected for motion, a linear trend and normalized to percentage signal change. The hemodynamic matrix H was based on the SPM Canonical HRF and the DCST was used to define the semiparametric matrix Φ. Here, the regularization parameter was selected to the universal threshold given by ˆ 2 log N λ σ = , where σ̂ is the MAD estimate